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Abstract It has long been recognized that aortic root elasticity helps to en¬ 
sure efficient aortic valve closure, but our understanding of the functional im¬ 
portance of the elasticity and geometry of the aortic root continues to evolve 
as increasingly detailed in vivo imaging data become available. Herein, we de¬ 
scribe fluid-structure interaction models of the aortic root, including the aortic 
valve leaflets, the sinuses of Valsalva, the aortic annulus, and the sinotubular 
junction, that employ a version of Peskin’s immersed boundary (IB) method 
with a finite element (FE) description of the structural elasticity. We develop 
both an idealized model of the root with three-fold symmetry of the aortic 
sinuses and valve leaflets, and a more realistic model that accounts for the 
differences in the sizes of the left, right, and noncoronary sinuses and corre¬ 
sponding valve cusps. As in earlier work, we use fiber-based models of the valve 
leaflets, but this study extends earlier IB models of the aortic root by employ¬ 
ing incompressible hyperelastic models of the mechanics of the sinuses and 
ascending aorta using a constitutive law fit to experimental data from human 
aortic root tissue. In vivo pressure loading is accounted for by a backwards 
displacement method that determines the unloaded configurations of the root 
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models. Our models yield realistic cardiac output at physiological pressures, 
with low transvalvular pressure differences during forward ffow, minimal re¬ 
gurgitation during valve closure, and realistic pressure loads when the valve 
is closed during diastole. Further, results from high-resolution computations 
demonstrate that IB models of the aortic valve are able to produce essentially 
grid-converged dynamics at practical grid spacings for the high-Reynolds num¬ 
ber ffows of the aortic root. 

Keywords aortic valve • ffuid-structure interaction • immersed boundary 
method • incompressible ffow • hyperelasticity • finite element method • finite 
difference method 


1 Introduction 

The aortic root consists of the three semilunar aortic valve cusps, the three 
sinuses of Valsalva, which are bulbous cavities positioned behind each valve 
leafiet, the aortic annulus, where the aortic valve leafiets attach to the aorta, 
and the sinotubular junction, where the sinuses merge into the ascending aorta. 
The healthy aortic valve acts to ensure unidirectional ffow of oxygenated blood 
from the heart to the tissues of the body, including to the heart itself via 
the coronary circulation. Diseases of the aortic valve can result in stenosis 
or regurgitation, and severe aortic valve disease is treated by repairing or 
by replacing the diseased valve with either a mechanical or a bioprosthetic 
valve [12]. Approximately 50,000 such procedures are performed each year 
in the United States [7,26,68]. Continuing advances in noninvasive in vivo 
imaging of blood ffow and tissue deformation, mechanical tests specifically 
targeted to biological tissues, and computer modeling and simulation enable 
integrative studies of the dynamics and mechanics of the aortic root [42] . Such 
work has the potential to improve both medical devices and surgical procedures 
used to treat patients with valvular heart disease, and also clinical approaches 
to patient risk assessment and treatment planning. 

Aortic compliance, and specifically the compliance of the aortic sinuses, 
has a fundamental role in the function of the aortic root. The sinuses act as 
reservoirs, storing blood during systole and then releasing it in diastole to fa¬ 
cilitate fiow in the coronaries [72]. The compliance of the sinuses also helps 
to ensure the proper closure of the aortic valve [4]. Further, the shape of the 
aortic root causes blood to circulate in the sinuses, generating vortices that act 
both on the valve leafiets and on the sinus walls. These vortices generate forces 
that facilitate effective valve closure and act as a regulatory mechanism that 
synchronizes closure [4]. Although this mechanism of efficient valve closure 
was first postulated by Leonardo da Vinci, there remains some debate in the 
surgical community about whether it is important to recreate the anatomic 
geometry of the native aortic root in aortic valve and root replacement proce¬ 
dures [13]. 

More recently, our understanding of the functional importance of the elas¬ 
ticity and anatomical geometry of the aortic root has evolved as increasingly 
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detailed data have been acquired on the in vivo deformations that occur during 
the cardiac cycle. For instance, it has been shown that the aortic root under¬ 
goes a multi-modal series of conformational changes even before the leaflets 
open [19]. It has been argued that these deformations act to reduce shear 
stresses on the valve leaflets and thereby prolong the life of the native valve 
leaflets [13]. Annular and commissural flexibility may be a key component 
in this interaction, and annular flexibility is lost with the implantation of a 
stented artificial valve, possibly reducing the lifetime of bioprosthetic valve 
leaflets. The annular commissures also participate in load sharing, which re¬ 
duces peak stresses on the aortic cusps immediately after valve closure. Devices 
such as the Medtronic 3f Aortic Bioprosthesis attempt to account for both the 
annular and the commissural flexibility of the native root, but despite the 
theoretical advantages offered by such designs, more traditional devices with 
rigid annuluses still dominate clinical practice. 

A challenge in modeling and simulating the mechanical response of arter¬ 
ies, including the aorta, is that the arterial walls are continuously subjected to 
substantial intraluminal pressures in their in vivo state. Thus, arterial geome¬ 
tries that are derived from in vivo imaging generally correspond to a loaded 
configuration. In addition, as first demonstrated by Fung [27] and by Vaish- 
nav and Vossoughi [73], even the unloaded configuration of the artery is not 
stress free. Instead, because of tissue growth and remodeling, the arterial wall 
includes residual stresses and is subject to axial tethering [11,28,42]. These 
factors all conspire to make determining the unloaded arterial geometry chal¬ 
lenging. Indeed, both the deformation due to intraluminal pressure loading 
and the residual stresses cannot be measured directly in vivo, and therefore 
these quantities must be estimated. Residual stresses may be determined from 
ex vivo experiments, whereas the zero-pressure configuration of blood vessels 
can be determined by solving the inverse elastostatic problem. 

Different approaches have been developed to solve the inverse elastostatic 
problem for complex arterial geometries derived from medical imaging data. 
Most of these strategies can be grouped in two broad categories. One group 
of approaches focuses on retrieving the initial deformation field and the initial 
stress field of the artery while keeping the imaging-derived geometry unaltered. 
An example of this type of approach is the modified updated Lagrangian for¬ 
mulation (MULF) introduced by Gee et al. [31], which is an incremental pre¬ 
stressing method that recovers the equilibrium configuration rather than the 
zero-pressure geometry. A second group of approaches determines an unloaded 
configuration of the vessel that, when subjected to intraluminal pressure load¬ 
ing, will inflate to match the imaging-derived geometry. Among the first to 
develop a stationary method to solve the inverse elastostatic problem were 
Govindjee and Mihalic [32]. This method was extended by Lu et al. [47] to 
arteries, but this approach does not appear to have been widely used in prac¬ 
tice, possibly because it requires a modification in the finite element (FE) solu¬ 
tion scheme that renders this method somewhat difficult to implement within 
commercial FE analysis software. In comparison, the backward incremental 
method developed by de Putter et al. [21] requires updating the deformation 
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gradient and thereby can be implemented with fewer modifications to the fi¬ 
nite element code. Perhaps the most straightforward approach to retrieving 
the zero-pressure configuration is the backward displacement method proposed 
by Bols et al. [6], which iteratively modifies the coordinates of the reference 
configuration rather than the deformation gradient tensor. We use the back¬ 
ward displacement method in this work to recover the unloaded configuration 
of an idealized model of the aortic root. 

While the methods to obtain an estimation of the zero-pressure configu¬ 
ration of a blood vessel can be applied to any geometry, determining residual 
stresses is a local, vessel-specific task that presently requires harvesting and 
then destroying the blood vessel. Computational approaches to accounting for 
residual stresses are usually based on an assumed ‘open’ configuration, which 
is subsequently ‘closed,’ thereby reversing the process that releases residual 
stresses [2,16]. This approach is difficult to apply to noncylindrical geome¬ 
tries and, in particular, cannot be readily applied to imaging-derived geome¬ 
tries [16]. A more generic approach to including residual stress in arteries is 
described in Alastrue et al. [2] and is based on a Kroner-Lee decomposition 
of the deformation gradient tensor. For simplicity, we do not consider residual 
stresses in the present work, although we do account for initial strains resulting 
from in vivo pressure loading. 

Most the methods surveyed above have been developed and applied to 
study aneurysmal arteries [31,47], although applications to normal vessels as 
well as other pathophysiological conditions are possible. Because arterial wall 
stress distributions are considered good predictors of aneurysm rupture, im¬ 
provements in the accuracy of computed wall stress distributions are expected 
to facilitate improvements in aneurysm rupture prediction. In contrast, al¬ 
though the important role of aortic root compliance in valve closure is well 
documented [17,18], many dynamic analyses of the aortic valve tend to focus 
on the mechanical behavior of the valve leaflets and often model the sinuses 
and the ascending aorta as rigid [35, 39, 75, 79] or as linear elastic materi¬ 
als [15,50], but experimental data [3] show that the sinuses and the ascending 
aorta exhibit a nonlinear hyperelastic behavior. Examples of nonlinear hyper¬ 
elastic constitutive models used in fluid-structure interaction simulations of 
the aortic root can be found in studies by De Hart et al. [20] and by Weinberg 
and Mofrad [76], which both adopt an arbitrary Lagrangian-Eulerian (ALE) 
approach. The application of ALE fluid-structure interaction schemes to mod¬ 
eling aortic valve dynamics is somewhat limited by the remeshing required by 
such methods as the structure deforms [69]. Attempts to tackle this limita¬ 
tion have included restricting the problem to a two-dimensional analysis [24] 
and the implementation of complex remeshing algorithms in three-dimensional 
analyses [14]. 

An alternative approach to modeling interactions between the blood and 
the aortic valve and root is offered by Beskin’s immersed boundary (IB) 
method [59], which was introduced to simulate cardiac valve dynamics [58] 
and was subsequently extended to model fluid-structure interaction in the 
heart and its valves [35-37,39,45,48,49,52-54,61]. This IB approach to fluid- 
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structure interaction is to describe structural stresses and deformations in La- 
grangian form, and to describe the momentum, viscosity, and incompressibil¬ 
ity of the coupled fluid-structure system in Eulerian form. Integral transforms 
with Dirac delta function kernels mediate interaction between Lagrangian and 
Eulerian variables. When discretized, the IB formulation of the equations of 
motion replaces these singular kernel functions with regularized kernels that 
are designed to ensure conservation of physical quantities such as force and 
torque when converting between Lagrangian and Eulerian forms. This numeri¬ 
cal approach allows the Lagrangian structural mesh to overlay the background 
Eulerian grid in an arbitrary manner, thereby avoiding the need to deploy dy¬ 
namic body-fitted grids. In addition, because the immersed structures move 
according to a common interpolated velocity field, the IB method offers an 
implicit contact model that prevents the valve leaflets from interpenetrating, 
even when subjected to substantial diastolic pressure loads [35,39]. Related 
sharp-interface IB methods have also been developed and used by Borazjani, 
Ge, and Sotiropoulos to simulate valvular dynamics for models of both rigid 
and flexible aortic valve prostheses [8,9]. Previous work demonstrated that 
three-dimensional IB models of the aortic valve can yield physiological cardiac 
output at realistic pressures [35,39]. However, we are aware of no previous 
study that demonstrates that the IB method yields reasonably well-resolved 
simulation results for flexible aortic valve models in three spatial dimensions 
at practical grid spacings. 

The primary aim of this study is to develop new fluid-structure interaction 
models of the aortic root that substantially extend earlier IB models of aortic 
valve dynamics [35,39] by including descriptions of the elasticity of the aor¬ 
tic sinuses and the ascending aorta. Herein, the aortic sinuses and proximal 
ascending aorta are modeled as an incompressible, isotropic, hyperelastic ma¬ 
terial with an exponential neo-Hookean strain-energy functional [22,42] that 
we fit to experimental data obtained by Azadani et al. [3] from human aortic 
root tissue samples. A separate finite element analysis is employed to esti¬ 
mate the unloaded aortic geometry using a method based on the backward 
displacement method of Bols et al. [6]. Our fluid-structure interaction simula¬ 
tions employ hybrid models of the aortic root that use a fiber-based description 
of the thin aortic valve leaflets, as done in previous studies [35,39], along with 
a finite element-based description [40] of the comparatively thick walls of the 
aortic sinuses and proximal ascending aorta. 

Two related aortic root models are considered in this work. Each captures 
the complex interactions between the flow, the thin flexible aortic valve cusps, 
and the deformable walls of the aortic sinuses and the ascending aorta. The 
first model employs a highly idealized anatomical geometry, in which the initial 
and reference configurations of the aortic sinuses and valve cusps are assumed 
to exhibit three-fold symmetry. The second model employs a more realistic 
description of the aortic root anatomy that accounts for the differences in the 
sizes of the left, right, and noncoronary sinuses. Both models are fully three 
dimensional in the sense that there are no symmetry conditions imposed on 
the subsequent fluid dynamics or structural kinematics. Indeed, because of 
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the relatively high Reynolds numbers of the systolic flow, symmetry breaking 
causes the three leaflets to move asynchronously even in the highly idealized 
aortic root geometry. 

Using the idealized model, we perform a refinement study to demonstrate 
that the present methods are able to yield essentially grid-converged results 
at practical grid spacings. Specifically, under grid refinement, only relatively 
small differences are observed in stroke volume, maximum flow velocity, and 
vessel distensibility and stress distributions. Relatively large differences remain 
in the details of the flow in the vicinity of the valve cusps, and in the details of 
the kinematics of the valve cusps. Higher spatial resolution is likely needed to 
resolve fully the fine-scale dynamics of the valve during systole. Nonetheless, 
a key contribution of this study is that it demonstrates, for the first time, that 
even bulk flow properties are reasonably resolved by IB models of aortic valve 
dynamics at practical spatial resolutions. 

We also compare the numerical predictions of the symmetric and asym¬ 
metric aortic root models. Previous studies have considered the influence of 
asymmetric leaflets, such as bileafiet aortic valves, in symmetric aortic root 
configurations [17,76]. In this work, we focus on assessing how symmetric and 
asymmetric geometries affect bulk flow hemodynamics and aortic wall mechan¬ 
ics. In particular, we aim to determine the extent that a symmetric geometry 
can be considered a reliable model of the aortic root, and to determine situ¬ 
ations where an asymmetric geometry gives a more physiologically accurate 
results. 

In our dynamic analyses, we pace the aortic root model to an essentially pe¬ 
riodic steady state using a prescribed, periodic left-ventricular pressure wave¬ 
form. A Windkessel model provides downstream loading for the aortic root. 
Notice that in these models, flow rates are not prescribed, but rather emerge 
from the computation. In all cases, physiological cardiac output is obtained 
at physiological driving and loading pressures, with low transvalvular pressure 
differences during forward flow, minimal regurgitation during valve closure, 
and realistic transvalvular pressure loads when the valve is closed during di¬ 
astole. 


2 Methods 

2.1 Immersed Boundary Formulation 

The immersed boundary formulation used herein describes fluid-structure sys¬ 
tems in which an elastic structure is immersed in a viscous incompressible 
fluid. It employs a Lagrangian description of the structural deformations and 
the stresses generated by those deformations, and an Eulerian description of 
the momentum, incompressibility, and viscosity of the coupled fluid-structure 
system. In the present model, we employ a fiber-based description of the thin 
valve leaflets, and we describe the aortic wall as an incompressible hyperelastic 
solid. Let q G 1/ C indicate Lagrangian curvilinear coordinates attached 
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to the valve leaflets, with q = (g, r), let X G ^ C indicate the Lagrangian 
material coordinate system of the aortic wall, with X = (X, F, Z), and let 
X G 17 indicate the Eulerian physical coordinates of the physical domain, with 
X = (x^y^z). The physical position of fiber point q at time t is given by 
0(^5 ^ the physical position of aortic wall material point X at time 

t is given by x(X,t) G 17. We use l7^(t) C 17 to indicate the (codimension 1) 
physical region occupied by the valve leaflets at time t, 17^ (t) C 17 to indicate 
the (codimension 0) physical region occupied by the solid-body model of the 
vessel wall, and l7^(t) = 17 \ 17^ (t) to indicate the physical region occupied by 
the fluid at time t. 

The equations of motion for the coupled fluid-structure system are: 


Du. . 

-Vp(x,t) + /uV^u(x,i) + f(x,t) + g(x,i), 

(1) 

V • u(x, t) = 

0 , 

( 2 ) 

f(x,t) = 

f F(q,t)(5(x-(/)(q,i))dq, 

Ju 

( 3 ) 

g(x,t) = 

[ Vx-P(X,t)(5(x-x(X,i))dX 

Jv 



- [ P(X,t)N(X)5(x-x(X,t))dA(X), 

( 4 ) 


JdV 


/ u(x,i)(5(x-</)(q,t))dx, 

JQ 

( 5 ) 


f u(x,i)(5(x-x(X,t))dx, 

Jn 

( 6 ) 


in which p is the mass density, p is the dynamic viscosity, u(x, t) is the Eulerian 
velocity field of the fluid-structure system, is the Eulerian pressure 

field of the fluid-structure system, f(x, t) is the Eulerian elastic force density 
generated by deformations to the fiber model of the valve leaflets, g(x, t) is 
the Eulerian elastic force density generated by deformations to the solid-body 
model of the vessel wall, F(q, t) is the Lagrangian elastic force density of the 
fiber model, P(X, t) is the first Piola-Kirchhoff elastic stress tensor of the solid 
model, N(X) is the outward unit normal along 9E, (5(x) = 6{x) 6{y) 6{z) is 
the three-dimensional Dirac delta function, and -b u • V(') is the 

material derivative. In the present work, we assume uniform mass density p 
and dynamic viscosity p. These assumptions are not essential, however, and 
versions of the IB method have been developed that permit the use of spatially 
varying mass densities [25,43,44,55,63,80] and viscosities [25,63]. 

Eqs. (1) and (2) are the Eulerian incompressible Navier-Stokes equations. 
Here, the momentum equation (1) is augmented by two Eulerian body force 
densities. The first of these, f(x, t), is determined by eq. (3) to be the Eulerian 
force density that is equivalent to the Lagrangian fiber force density F(q, t). 
The second body force in (1), g(x,t), is determined by eq. (4) to be the 
Eulerian force density equivalent to the Lagrangian description of the forces 
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generated by the solid-body model of the vessel wall, which are expressed in 
terms of P(X,t). Notice that f(x, t) is a singular force density supported on 
By contrast, g(x, t) is supported on Away from g(x, t) is 

nonsingular, although g(x, t) is singular along df2^{t) if PN 7 ^ 0 ; see eq. (4). 

Eqs. (5) and ( 6 ) determine the motions of the immersed structures from 
the Eulerian material velocity field u(x, t). Because of the presence of viscosity 
throughout the fluid-solid system, u(x, t) is continuous, and thus (5) and ( 6 ) 
are equivalent to 


— (q,i) = u(</)(q,t),i), (7) 

= (8) 

Because V • u(x, t) = 0, the immersed solid is automatically treated as incom¬ 
pressible in this formulation. Specifically, at least in the continuum formula¬ 
tion, there is no need to include terms in the elastic strain-energy functional 
to penalize compressible deformations. 

In practice, we use a standard finite element method to approximate 
the deformations and stresses of the vessel wall. To do so, it is useful to adopt 
a weak formulation for the structural equations, so that 


g(x,t)= [ G(X,i)<5(x-x(X,i))dX, 

Jv 

[ G(X,t) ■V(X)dX = - [ P(X,i) : VxV(X)dX, W(X), 
Jv Jv 


^(X,t)=U(X,t), 


[ U(X,t)-V(X)dX= [ u(x(X,i),t)-V(X)dX, 
Jv Jv 


( 9 ) 

( 10 ) 

( 11 ) 

( 12 ) 


in which G(X,t) is the Lagrangian elastic force density of the aortic wall, 
U(X, t) is the Lagrangian velocity field of the aortic wall, and V(X) is an ar¬ 
bitrary Lagrangian test function that is not assumed to vanish on dV. In the 
continuum setting, eqs. (4) and (9)-(10) are equivalent definitions for g(x, t); 
however, these formulations lead to different numerical schemes when dis¬ 
cretized, and only eqs. (9)-(10) lead directly to a standard nodal finite element 
structural discretization. Eurther, in the continuum equations, eqs. (11)-(12) 
are equivalent to eq. ( 6 ) and to eq. ( 8 ). 


2.2 Valve Leaflets 
2.2.1 Geometry 

As in earlier work [35,39], the leaflet geometry is determined by the mathe¬ 
matical theory of the fiber architecture of aortic valve leaflets developed by 
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Peskin and McQueen [60] . This theoretical model of valve geometry and fiber 
architecture derives the shape of the leaflets from their function, which is to 
support a uniform pressure load during diastole, when the valve is closed. Each 
valve leaflet is composed of two families of fibers that are orthogonal to each 
other. One family of fibers runs from commissure to commissure, and the other 
runs from the bottom scalloped edge of the aortic sinus to the free edge of the 
valve leaflet. The leaflets are defined in a closed and loaded configuration; see 
fig. 1. In this configuration, the radius of each leaflet, measured from the tip of 
the valve to the end of the belly region, is 1.45 cm, and the coapting portion of 
each leaflet is 0.97 cm tall. The scale of the leaflets is based on measurements 
of human aortic roots [62,71]. The height of the coapting portion of the leaflet 
is chosen to be similar to that of the real valve while enabling the model valve 
to support a physiological pressure load when closed. 


2.2.2 Mechanical response 


The mechanical behavior of the leaflets is described by a strain-energy func¬ 
tional E = E[0(-,t)] described previously [35,39]. Briefly, let the curvilinear 
coordinates q = (^, r) be chosen so that a fixed value of q labels an individual 
fiber, and so that r runs along the fibers. Consequently, {d(p/dr)/\\dcj)/dr\\ 
is the unit fiber tangent vector. The total elastic energy E is the sum of a 
stretching energy E^ and a bending energy Eb, 


E — Eg ^ -Eb, 



dq, 

dr‘^ 


2 

dq. 


(13) 

(14) 

(15) 


Eq. (14) accounts for the total stretching energy associated with the fibers, in 
which Ss is a spatially inhomogeneous local stretching energy with a quadratic 
length-tension relationship, and eq. (15) accounts for the total bending energy 
associated with the fibers, in which Cb is a spatially inhomogeneous bending 
stiffness and 0 is the reference configuration, which is taken to be the initial 
configuration. Because the valve leaflets are modeled as thin elastic surfaces, 
the bending resistant energy allows the model to account for the thickness of 
the real valve leaflets. Larger values of Cb model a thick, stenotic valve, whereas 
smaller values of Cb model a thin, flexible valve. The resulting Lagrangian body 
force density is the Erechet derivative of E. 

Valve leaflet parameters are empirically determined for a fixed leaflet dis¬ 
cretization by choosing the leaflet stiffness so that the leaflet supports a dias¬ 
tolic load of approximately 80 mmHg with 10% strain in the family of fibers 
running from commissure to commissure. This leads to a commissural fiber 
stiffness of 7.5e6 dyne/cm. We further assume that these commissural fibers 
are a factor of 10 stiffer than the radial fibers that run orthogonal to the com¬ 
missural fibers [65]. The bending stiffness is chosen to be approximately the 
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Fig. 1: Geometrically idealized model of the aortic root. We use Peskin and 
McQueen’s theoretical model of the collagen fiber architecture of aortic valve 
leaflets [60] to determine the shape of the leaflets for (a) a symmetric and (b) an 
asymmetric aortic valve. We use measurements of human aortic roots [62,71] 
to determine the geometry of the outflow tract, the sinuses, and the ascending 
aorta. Panel (c) shows a side view and a top view of the symmetric aortic 
root model, and panel (d) shows a side view and a top view of the asymmetric 
aortic root model. 


smallest value that yields valve leaflets that successfully coapt at the end of 
systole. These material parameters are chosen only to yield a functional valve, 
and are not represented as corresponding to those of an actual valve cusp. 
Nonetheless, this model of the valve leaflets does account for key structures 
of actual valve cusps, including the prominent fibers that run from commis¬ 
sure to commissure and the 10-to-l anisotropy of real valve leaflets [65]. See 
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refs. 39 and 35 for further details. In preliminary studies, a sensitivity analysis 
was conducted for the bending stiffness of the aortic leaflets, and it was found 
that variations on the order of ± 25% in the bending stiffness have little effect 
on the physiological predictions of the models (not shown). 


2.3 Aortic Wall 
2.3.1 Geometry 

An idealized model of the aortic root is constructed, as in previous work [35,39]. 
The dimensions of this model are based on measurements by Swanson and 
Clark [71] of human aortic roots collected at autopsy, and the geometry of 
the model is based on measurements by Reul et al. [62] from angiograms. The 
diameter of the aortic portion of the model is 3 cm, whereas the diameter of 
the sinus region, measured from a commissura to the center of a sinus, is 3.5 
cm. The overall length of the model is of 10 cm, and the distance between 
the annulus and the aortic flow outlet is 7.75 cm. The thickness is constant 
throughout the model and is 2 mm [67]. In our simulations, we employ a 
semi-rigid model of the left-ventricular outflow tract while allowing for a fully 
flexible description of the aortic sinuses and the ascending aorta. The aortic 
annulus, which is the scalloped line of attachment between the valve leaflets 
and the aortic sinuses, was considered to be the lower boundary of the sinuses; 
see fig. 1. The aortic model used in this study differs from previous work 
[35, 39] in that here we treat the aortic sinuses and ascending aorta as an 
incompressible hyperelastic material, whereas in earlier studies it was treated 
as a rigid structure. 

In our simulations, we consider both symmetric and asymmetric models 
of the aortic root. In the symmetric model, we assume that the initial and 
reference configurations of the aortic root exhibit three-fold symmetry, so that 
each sinus and the corresponding leaflet occupies an angle of 120°. For the 
asymmetric model, the symmetric geometry is modified so that each sinus, 
along with the corresponding leaflet, covers a different angle. Based on the 
work of Berdajs [5], the average angle covered by each leaflet was found to 
be 136.22° for the right coronary sinus, 122.48° for the noncoronary sinus, 
and 101.3° for the left coronary sinus. Other characteristic dimensions of the 
leaflets, such as the radius or the height of the coapting region, were taken to 
be the same as the symmetric model; see fig. 1. 

Because the geometry of the valve leaflets are defined in a closed and 
loaded configuration, it is also necessary to specify the geometry of the vessel 
wall in a loaded configuration, and to compute the corresponding unloaded 
configuration. This configuration will depend on both the prescribed loaded 
configuration and the constitutive model. To determine the unloaded configu¬ 
ration, we employ an iterative method described by Bols et al. [6] detailed in 
sec. 2.5.4. 
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2.3.2 Mechanical response 

We describe the elasticity of the aortic sinuses and ascending aorta using a 
hyperelastic constitutive model fit to biaxial tensile test data collected by 
Azadani et ah [3] from tissue samples from human aortic sinuses. Because the 
data reported by Azadani et al. indicate a nearly isotropic material response 
(see fig. 2 ), in this work we model the material response of the aortic sinuses 
and ascending aorta using an isotropic strain-energy functional W with an 
exponential stress-strain relation, 

^ [exp( 6 (/i - 3)) - 1], (16) 

in which Ii = /i(C) = tr(C) is the first invariant of the right Cauchy-Green 
strain tensor C = F^F, and F = 9x/9X is the deformation gradient tensor 
associated with the deformation mapping x • Material parame¬ 

ters were obtained by least-squares fits to experimental data using MATLAB 
(Mathworks, Inc., Natick, MA, USA). The first Piola-Kirchhoff stress tensor 
F is determined from W via 

P = ^ - Ps F-^ = cexp( 6 (/i - 3)) F - ^ F”^, (17) 

in which ps is a structural pressure-like term that is chosen to improve the 
accuracy of the stress predictions of the method [29]. As in earlier work [29], 
we compute Ps via 


Ps = cexp(&(Ji - 3)) - /3s log(/ 3 ), (18) 

in which /s = / 3 (C) = det(C) is the third invariant of C and (3s = 2.5e4 kPa. 
Thus, F = 0 for F = I, and the structural model provides an energetic penalty 
for any compressible deformations. 


2.4 Driving and Loading Conditions 

A left-ventricular pressure waveform adapted from human clinical data of 
Murgo et al. [56] is prescribed at the upstream inlet of the left-ventricular 
outflow tract to drive flow through the model aortic root, and downstream 
loading conditions are provided by a three-element Windkessel model by Ster- 
giopolus et al. [70] fit to the clinical data of Murgo et al. [56]. As in earlier 
work [35,39], zero pressure boundary conditions are imposed on the boundary 
of the fluid-filled region exterior to the aortic root model. Coupling between 
the detailed IB model and the reduced circulation models is also done as in 
earlier work [35,39]. 
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2.5 Numerical Methods 

We use a locally-refined staggered-grid discretization of the Eulerian equa¬ 
tions along with a finite difference-based discretization of the fiber model of 
the valve leaflets and a finite element-based description of the continuum ves¬ 
sel wall model. Our approach is similar to the spatially adaptive IB scheme 
described previously [35], except that here we also employ a finite element- 
based description of the aortic wall, using an approach introduced by Griffith 
and Luo [40]. 


2.5.1 Eulerian and Lagrangian spatial discretizations 

The physical domain i? is discretized using a locally-refined Cartesian grid, but 
for simplicity we describe only the uniform-grid version of this method; details 
on the adaptively refined discretization are provided by ref. 35. Let (i,j,/c) 
index the Cartesian grid cells, and let = ((^ + (j + 

indicate the position of the center of grid cell {i^j, k), in which h is the Carte¬ 
sian grid meshwidth and Z\x = is the Cartesian grid cell volume. The 
Eulerian velocity field u = n, w) is approximated at the center of each face 
of the Cartesian grid cells in terms of the velocity component that is normal 
to that face: u is approximated at the locations x^_ i v is approximated at 

the locations x^ ^_i and w is approximated at x^ ^- . The Eulerian force 

densities f and g are approximated in the same staggered-grid fashion. The 
Eulerian pressure p is approximated at the centers of the Cartesian grid cells. 

The deformations and forces associated with the valve leaflets are approxi¬ 
mated on a fiber-aligned curvilinear mesh. Let I index the nodes of this mesh, 
let (f)i and F/ indicate the current position and Lagrangian force density of 
node /, respectively, and let Z\q^ indicate the area fraction (quadrature weight) 
associated with node 1. F/ is computed from using a finite difference ap¬ 
proximation to the fiber force density; see refs. 39 and 35 for details, including 
relevant finite difference formulae. 

The deformations, stresses, and resultant forces associated with the aortic 
wall are approximated on a volumetric Lagrangian mesh. Let e index the 
elements of this mesh, with W indicating the volume associated with element 
e, so that V = LJeW, let m index the mesh nodes, and let indicate the 

interpolatory Lagrangian finite element basis function associated with node 
m. The structural deformation and elastic force density are approximated in 
a standard manner via 


= (19) 

m 

G(X,t) = ^G„{t)^„(X), (20) 


in which are the nodal positions and G^(t) are the nodal force densi¬ 

ties. The deformation gradient F is computed by directly differentiating the 
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approximation to x(X,t), and P is evaluated from the approximation to F. 
The nodal values G^(t) are determined so that 

[ G(X,t)(^,(X)dX = - [ F(X,t)Vx^n(X)dX (21) 

Jv Jv 

for all n. This leads to a linear system of equations that defines the nodal force 
densities in terms of the nodal deformations. In practice, these integrals are 
evaluated element-by-element using Gaussian quadrature rules that are exact 
for the left-hand side but are generally only approximate for the right-hand 
side. 


2.5.2 Lagrangian-Eulerian interaction 

To couple the Lagrangian and Eulerian discretizations, we employ approxi¬ 
mations to the integral transforms of the continuum equations that replace 
the singular Dirac delta function (5(x) = 5{x) 5{y) 5{z) by a regularized delta 
function 5h(x) = 5h(x) 5h(y) 5h(z). In our computations, we construct the 
three-dimensional regularized delta function either by using the four-point 
one-dimensional regularized delta function of Peskin [59], or by using a broad¬ 
ened version of this function that has a spatial extent of 8 meshwidths. The 
same regularized delta function is used for both the thin model of the valve 
leaflets and the volumetric model of the vessel wall. 

For the leaflet model, we employ a coupling approach frequently used with 
the IB method. The Lagrangian forces F = associated with the 

leaflets are converted into equivalent Eulerian forces f = (/^, /^, /^) via 


1 

(22) 

1 

(23) 

= E - <t>i)Ml- 

(24) 


i 


This amounts to using the trapezoidal rule to discretize the integral transform 
(3). We employ the compact notation 

f = 5'[0]F (25) 


to express this force-spreading operation. The Eulerian velocity u = {u^v^w) 
is used to determine the dynamics of the physical positions (j) = (0^,0^,^^) 
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of the material points of the valve leaflets via 


- (t>l) (26) 

- 4>i) ^x, (27) 

= J2 - <t>l) (28) 

We use the notation 

^(f) = 'R}[cl)]u (29) 


to express this velocity-restriction operation. Notice that by construction, the 
force-spreading and velocity-interpolation operators are adjoints, i.e., = 

The Lagrangian-Eulerian coupling scheme used for the vessel model is 
based on an approach developed by Griffith and Luo [40]. This approach is 
first to construct a force-spreading operator, and then to determine a velocity- 
restriction operator that is the adjoint of the force-spreading operator. Briefly, 
for each element e, a set of quadrature points Q{e) is determined, and for 
each quadrature point q G Q(e), X® indicates the reference coordinates of 
the quadrature point, and uj^ indicates the quadrature weight associated. The 
Lagrangian force density G = {G ^is converted into the equivalent 
Eulerian force density g = by discretizing eq. (9), with (5(x) re¬ 

placed by Sh(x), using this quadrature rule, i.e. 


X 

dik 

^id-^,k 


Gyx^,^)<5;,(x,_l,^.,-x(x^,^)):^,^ 

(30) 

e geQie) 



(31) 

6 q^Q{e) 



(32) 


e qeQ{e) 


Notice that these formulae take advantage of the fact that we may evaluate 
the approximations to G and x arbitrary Lagrangian locations. Specifically, 
we are not restricted to evaluating G and x f^he nodes of the finite element 
mesh. In practice, the quadrature rules are dynamically determined to ensure 
that the physical spacing of the quadrature points is on average one-half of a 
Cartesian meshwidth, even under the presence of extremely large structural 
deformations. We use the notation 




( 33 ) 


to denote the force-spreading operator. 
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To determine the corresponding velocity-restriction operator, we introduce 
a Lagrangian velocity field U = ([/, V, IT), which is defined in terms of nodal 
velocities via 

U(X,t) = ^U„{t)^„(X). (34) 

m 

This velocity field is required to satisfy, for all n, 

e qeQ(e) 




e qeQ(e) \i,j,k 


(35) 


^ ^ y(x^,t)^„(xp 

e qeQ{e) 




= T] T] - x(^‘q,t))Ax (^„(X®)w®, 

e qeQie) \i,j,k ) 

(36) 


e qeQ{e) 




= Z] H lj2'^iJ,k-i^h{^iJ,k-^-X(^q,t))Ax\ipniXl)uj^g. 
e qeQie) \i,j,k ) 

(37) 


This yields a system of linear equations to be solved for the nodal values of U. 
Notice that eqs. (35)-(37) correspond to a discrete component-wise approxi¬ 
mation to eq. (12). We use the notation 

^X = U = 7^-Mu. (38) 

Notice that 7^^[x] is a nonlocal operator in the sense that it requires the 
solution of a system of linear equations. Although not shown here, it is the 
case that See Griffith and Luo [40] for further discussion. 


2.5.3 Time stepping 

The time stepping scheme is similar to that of Griffith and Luo [40]. Let At 
indicate the (uniform) time step size, and let [T,T+^] = [nZ\t, (n + l)4\t] 
indicate the time interval. In this subsection, superscript indices always 
indicate the time step number. The state variables u, 0, and % are defined at 
integer time steps. The pressure p, which is not a state variable of the system, 
is defined at half-integer time steps. Approximations to the state variables 
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and to derived quantities such as the Lagrangian forces are also evaluated at 
half-integer time steps. 

The time stepping scheme proceeds as follows. First, predicted intermedi¬ 
ate approximations to the structural deformations 0 and x time are 
determined via 


At/2 

_ X^ 

At/2 


[(/)”] u", 


(39) 

(40) 


This is an application of the forward Euler time stepping scheme. Lagrangian 
force densities and are determined from the predicted structure 

configurations, and these Lagrangian force densities are spread to the Carte¬ 
sian grid via 


=5'[(/)”+^]F”+b (41) 

=5"'[x"+5]G"+i (42) 


Next, we solve the discretized incompressible Navier-Stokes equations for 
and p ’^+2 via a Crank-Nicolson Adams-Bashforth scheme. 


P 


/u^+i 

V M 


[u ■ = -VhP^+i + pVl 


+ f"+5 +g«+i 


Vft • u"+i = 0, 


(43) 

(44) 


with 

[u ■ (45) 

In the uniform grid case, the discrete operators Vh * , and are stan¬ 
dard second-order accurate staggered-grid finite difference approximations to 
the gradient, divergence, and Laplace operators, respectively [34]. We com¬ 
pute u • V/iU via a staggered-grid version of the piecewise parabolic method 
(PPM) [34]. On locally-refined Cartesian grids, these discrete operators are 
straightforward extensions of the uniform grid operators [35]. Solving eqs. (43) 
and (44) for and requires the solution of the time-dependent 

Stokes equations. The Stokes equations are solved using an iterative Krylov 
method with a multigrid preconditioner derived from the classical projection 
method [10,34]. 

Finally, we determine approximations to the structural deformations at the 
end of the time step via 


0 ^+^ - 0 ^ 

~At 

^n+l _ 


72.‘[(^"+5]u”+t 

(46) 


(47) 


At 
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with 


u 


n+i 


U^+1 + 
2 


(48) 


This is an application of the explicit midpoint rule to the structural dynamics. 

In the initial time step, we do not have a value for as required to 

evaluate the Adams-Bashforth approximation to the convective term. Conse¬ 
quently, in the initial time step, we employ a two-step Runge-Kutta scheme. 

Within each time step, coupling between the detailed fluid-structure in¬ 
teraction model and the reduced Windkessel model of the circulation is done 
as described previously [35,39]. In this approach, we perform only a single 
Stokes solve per time step, except in the initial time step, when we perform 
two Stokes solves within the context of a two-step Runge-Kutta scheme. 


2 . 5.4 Backward displacement method 


To determine the unloaded configuration of the vessel wall, we use an itera¬ 
tive backward displacement method described by Bols et al. [6] implemented 
using a custom MATLAB script that interfaces with the ABAQUS (Dassault 
Systemes Simulia Corp., Providence, RI, USA) finite element analysis soft¬ 
ware. Because the constitutive model (16) is not provided by ABAQUS, this 
model was implemented using a UHYPER subroutine. 

Let denote positions of the finite element mesh nodes in the prescribed 
loaded configuration, let denote the corresponding nodal positions in the 
computed reference configuration after i iterations of the backward displace¬ 
ment algorithm, and let x\n — denote nodal positions in the computed 

loaded configuration when X* is used as the unloaded reference configuration. 
Notice that x^, X^, and xln defined at the nodes of the finite element 

mesh. 

The algorithm is straightforward: First, we initialize X^ := x (i.e., we use 
the deformed coordinates as an initial guess for the unloaded configuration). 
Then, given X\ X*^^ is determined by 

X*+1 ;=x' + (x- V)- (49) 


Notice that if we define the displacement from the computed reference config¬ 
uration to the deformed coordinates at step i by D* = %* — X*, then (49) is 
equivalent to 

:= X - D\ (50) 

Thus, the reference coordinates are determined via a backward displacement 
from the deformed configuration. If this process converges, ^ x and we 
obtain reference coordinates X such that x[X] = x. 

The convergence of this iterative process is assessed in terms of the discrete 
norm of the residual x — i.e. 


r* = ||x-x[Xi||2 



(51) 
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The routine was judged to have converged once the residual was less than 
0.1% of the average vessel radius and the change in the residual between two 
consecutive iterations was less than 0.015% of the average radius. In general, 
obtaining convergence may require the use of damping, or of gradually increas¬ 
ing the loading pressure; however, in practice, we have found this algorithm 
to be robust, and neither damping nor incremental loading were required for 
the analyses performed herein. 

2.5.5 Discretization 

In our dynamic fluid-structure interaction simulations, the computational do¬ 
main is taken to be a 10 cm x 10 cm x 10 cm cubic region discretized us¬ 
ing a three-level locally refined grid. For the purposes of a mesh convergence 
study, two different fine grid spacings are used: a coarser one corresponding to 
h = 0.78125 mm, which is equivalent to a uniform 128 x 128 x 128 Eulerian 
discretization, and a finer one corresponding to h = 0.390625 mm, which is 
equivalent to a uniform 256 x 256 x 256 Eulerian discretization. The curvi¬ 
linear mesh used to discretize the valve leaflets has a physical grid spacing 
of approximately h/4 for the coarser Eulerian discretization, and to h/2 for 
the finer discretization. The volumetric mesh used to discretize the vessel wall 
uses 23,400 8-node (trilinear) hexahedral elements for both the symmetric and 
the asymmetric model. The average mesh aspect ratio was 1.93 ± 0.25 for the 
symmetric root model and 1.94 ± 0.34 for the asymmetric one. A preliminary 
convergence study in ABAQUS was performed to ensure mesh convergence of 
the structural model for both the backward displacement method and the EE 
analysis. We use the standard four-point delta function of Peskin [59] for the 
coarser simulations, and a broadened 8-point version of this kernel function for 
the higher resolution simulations, so that the physical extent of the regular¬ 
ized delta function is the same for both spatial resolutions.^ Consequently, the 
valve leaflets have the same effective thickness in for both spatial resolutions. 
A uniform time step size of At = 9.94369e-6 s was used with the coarser Eu¬ 
lerian discretization, whereas a uniform time step size of At = 1.40625e-5 was 
used for the finer Eulerian discretization. These time step sizes are empirically 
determined to be within a factor of ^/2 of the largest time step that satisfies 
the stability restriction of our explicit time integration method. 


3 Results 

3.1 Constitutive Parameters 

The constitutive model (16) used to describe the mechanical response of the 
aortic sinuses and ascending aorta is fit to data collected biaxial tensile tests 
of human aortic root tissue samples reported by Azadani et al. [3], yielding 

^ Similar convergence studies, although for substantially different applications, are de¬ 
scribed in refs. 38 and 29. 
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Stretch 


Fig. 2: Biaxial tensile test data obtained from human aortic root tissue samples 
by Azadani et al. [3] and our fit of the isotropic hyperelastic constitutive model 
(16) to these data. 


model parameter values of c = 12.8 kPa and b = 6.9; see fig. 2. Recall that 
these material parameters are used in the aortic sinuses and ascending aorta, 
whereas the outflow tract is treated as an essentially rigid structure. 


3.2 Unloaded Geometries 

The unloaded configuration is determined by assuming that the loaded con¬ 
figuration corresponds to a pressure load of 80 mmHg. In fig. 3, panels (a) 
and (b) show the prescribed loaded geometry, and panels (c) and (d) show 
the computed zero-pressure configuration. Fig. 4 shows the convergence his¬ 
tory of the backward displacement method. When the computed zero-pressure 
geometry is subjected to a 80 mmHg pressure load, the deformed geometry 
agrees with the prescribed geometry to within approximately 15 /im; see fig. 5. 
The discrepancy between the prescribed loaded geometry and the computed 
loaded geometry is greatest at the commissures and smallest along the straight 
portion of the vessel. Similar results are obtained for both the symmetric and 
asymmetric models. 






residual (m) 
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(a) (b) (c) (d) 


Fig. 3: Views of the (a) loaded and (b) unloaded symmetric geometry, and of 
the (c) loaded and (d) unloaded asymmetric geometry. 



Fig. 4: Convergence history of the inverse displacement procedure. 
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f 


(f) 

Fig. 5: Contour plots showing the norm of the distance between the prescribed 
and computed loaded geometries. Panels (a) and (b) show the symmetric ge¬ 
ometry, and panels (c)-(f) show the asymmetric geometry, with panels (c) and 
(d) showing the right sinus and panels (e) and (f) showing the left and non¬ 
coronary sinuses. Recall that the unloaded geometry is determined from the 
loaded geometry. The computed loaded geometry is the result of applying a 
prescribed pressure load to the computed unloaded geometry. 
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3.3 Fluid-Structure Interaction Analyses 
3.3.1 Symmetric model 

A dynamic fluid-structure interaction analysis is performed that includes four 
cardiac cycles, including an initialization cycle and three cycles in which the 
model has attained periodic steady state. (Data are reported only for the 
beats in which the model has reached steady state.) Fig. 6 shows the pre¬ 
scribed left-ventricular driving pressure along with the computed aortic load¬ 
ing pressure determined by the coupled Windkessel model and the computed 
flow rate through the symmetric aortic root. We emphasize that the flow rate 
is not imposed in the model, but rather is predicted by the fluid-structure 
interaction analysis. For the simulations using the coarser Cartesian spac¬ 
ing, h = 0.78125 mm, stroke volume is 96.2 ml, peak flow rate is 591.5 
ml/s, and cardiac output is 6.5 1/min. For the finer Cartesian grid spacing, 
h = 0.3906 mm, stroke volume is 100.2 ml, peak flow rate is 643.8 ml/s, and 
cardiac output is 6.8 1/min. Both sets of simulation results are within 15% 
of the hemodynamic parameters of the subject-specific clinical data used to 
determine the driving pressure and the Windkessel model parameters [56,70]. 
Specifically, the driving pressure and the Windkessel model parameters used 
in this work were based on clinical measurements from a subject with stroke 
volume of 100 ml, peak flow rate of 560 ml/s, and cardiac output of 6.8 1/min, 
as reported by Murgo et al. [56]. The simulations also capture the incisura in 
the flow profile corresponding to the backward flow generated by the leaflets. 
The backward flow volume is 1.77 ± 0.11 ml, or 1.8% of the stroke volume, for 
the coarser grid spacing, and 1.52 ± 0.059 ml, or 1.5% of the stroke volume, 
for the finer grid spacing. The net forward flow volume (i.e., into the aorta) 
following leaflet coaptation is 0.02 ± 0.01 ml. This forward flow results from 
interactions between the elastic aortic root and the dynamic downstream load¬ 
ing pressure. As the downstream pressure falls during diastole, the load on the 
aortic wall and leaflets decreases and allows the aortic root to gradually push 
fluid into the distal aorta. 

The model captures the opening and closing dynamics of the valve at both 
Cartesian grid spacings. Fig. 7 shows the leaflets being pushed apart along 
with the formation of the vortices that help the leaflets to close at the end 
of systole. Fig. 8 shows that although there is some net flow towards the left 
ventricle as the valve closes, once the valve closes, no further flow is allowed 
between the aorta and the left ventricle. 

Fig. 9 shows the distribution of the structural displacements in the aortic 
root, and fig. 10 shows displacement as a function of time at four selected 
locations: the center of the sinus; the sinotubular junction; the commissure; and 
the ascending aorta. Notice that after the first cardiac cycle, the displacements 
are essentially periodic in time. Cartesian grid resolution has a minimal effect 
on the range of displacements except for the displacements measured at the 
commissure. The largest amplitude is observed in the ascending aorta (0.58 ± 
0.05 mm) and at the sinotubular junction (0.52T0.02 mm). A measure of aortic 
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(a) 



(b) 

Fig. 6: Results of a dynamic fluid-structure interaction analysis of the sym¬ 
metric aortic root over the final three simulated cardiac cycles using a coarser 
Cartesian grid (effective grid resolution of = 128) and a finer Cartesian grid 
(effective grid resolution of N = 256). (a) Prescribed left-ventricular driving 
pressure and computed aortic loading pressure, and (b) computed flow rate, 
in which the dashed vertical lines highlight the moment that the valve leaflets 
open and the solid vertical lines highlight the moment that the valve leaflets 
close. Large oscillations indicate the reverberations of the aortic valve leaflets 
upon closure. 
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Fig. 7: Valve leaflet configurations during valve opening shown at equally 
spaced times during the fourth simulated cardiac cycle for (a) the coarser 
Cartesian grid (effective grid resolution of N = 128) and (b) the finer Carte¬ 
sian grid (effective grid resolution of N = 256). Panels (c) and (d) show the 
corresponding axial velocity profiles. Panel (e) shows the computed flow rates, 
with the dashed lines indicating the time instants shown in panels (a)-(d). 


stiffness is the radial pulsation, i.e., the ratio of the difference between systolic 
radius and diastolic radius to the average radius [57], which in the present 
analysis was approximately 3.9% for the coarser Cartesian grid spacing and 
2.9% for the finer Cartesian grid spacing. These values are similar to the value 
of 2.5% reported in refs. 42 and 57 for the ascending aorta. The distribution 
of the displacements was also in good agreement with literature data, with 
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(e) 


Fig. 8: Similar to fig. 7, but here showing the dynamics of the closure of the 
aortic valve during the fourth cycle. In this case, the solid lines indicate the 
instants shown in the figure panels. 


the displacements in the aortic root greater in the sinotubular junction, and 
at the commissures than in the sinuses, as reported also by Lansac [46]. 

The distribution of the maximum principal stress in the aortic root follows 
a distribution similar to that of the displacement, with the maximum stresses 
concentrated in the sinotubular junction, just above the sinuses, as shown in 
fig. 11. Stress in the aortic root follows the pulsatile waveform of the pressure; 
see fig. 12. Both Cartesian grid spacings yield similar stress distributions; 
compare fig. 11 panels (a) and (b), but fig. 12 shows some differences along 
the sinotubular junction and within the sinuses. The maximum principal stress 
in the sinuses is 20% of that observed in the rest of the aortic root. The 
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(c) 



Displacements (mm) 



Fig. 9: Displacement contours along the aortic root model. Panel (a) shows 
displacement contours along the axial section of the aortic wall during valve 
opening, shown as the same time instants as in fig. 7, for the coarser Cartesian 
grid spacing. Panel (b) is similar, but here shows results for the finer Cartesian 
grid spacing. Panels (c) and (d) are similar, but here show the displacement 
patterns during valve closure. 


distribution, and the range of the maximum principal stress predicted by our 
model, is in agreement with other analyses available in literature, in particular 
with the 220 kPa peak of the maximum principal stress in the sinotubular 
junction reported by Conti et al. [15]. 
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(a) (b) 



(c) (d) 

Fig. 10: Displacements of the deformable aortic root model at four locations: 
(a) the straight portion of the aorta; (b) just above one of the commissures; 
(c) the sinotubular junction; and (d) the middle of the sinus. 


3.3.2 Asymmetric model 

Similar fluid-structure interaction simulations are performed using the asym¬ 
metric model using only the coarser Cartesian grid spacing oih = 0.78125 mm. 
Fig. 13 shows the prescribed left ventricular driving pressure and the computed 
aortic pressure along with the computed flow rate generated by the model. 
Peak flow rate in this model is 581 ml/s, which is reduced compared to the 
symmetric studies and is within 5% of the subject-specific peak flow rate of 
560 ml/s [56]. Stroke volume and cardiac output are 95.1 ml and 6.46 1/min, 
respectively, which are both within 6% of the respective clinical values of 100 
ml and 6.8 1/m [56]. 

The details of the valve kinematics are also different with the asymmetric 
model, but the model is still able to capture the leaflets being pushed back, 
the flow inversion, the vortices in the sinuses, and the valve closure, which 
does not allow regurgitation during diastole; see figs. 14. 

Overall, the distribution of the displacements appeared similar to the sym¬ 
metric case, with the ascending aorta and the sinotubular junction moving 
more than the rest of the aortic root, as shown in fig. 15; however, tracking 
the displacements over time at specific locations for the right, left, and non¬ 
coronary sinus shows that the left sinus moves 20% less than the right sinus 



































Immersed boundary-finite element model of FSI in the aortic root 


29 



Max Principal Stress (kPa) 



Fig. 11: Maximum principal stress contours along the aortic root model. Panel 
(a) shows displacement contours along the axial section of the aortic wall 
during valve opening, shown as the same time instants as in fig. 7, for the 
coarser Cartesian grid spacing. Panel (b) is similar, but here shows results for 
the finer Cartesian grid spacing. Panels (c) and (d) are similar, but here show 
the displacement patterns durning valve closure. 


and 10% less than the noncoronary sinus. In addition, fig. 15 shows that the 
left sector of the aorta maintains a restricted motion compared to the right 
and noncoronary sectors of the aorta. This behavior can be observed also in 
the aorta in which the radial pulsation ranges from 1.9% for the left sector of 
the aorta, to 3.9% for the right sector of the aorta. The average displacement 
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(a) (b) 




(c) (d) 

Fig. 12: Maximum principal stress of the deformable aortic root model at 
four locations: (a) the straight portion of the aorta; (b) just above one of the 
commissures; (c) the sinotubular junction; and (d) the middle of the sinus. 


amplitude in the sinotubular junction is 0.49 mm for the right and noncoro¬ 
nary sinus and 0.36 mm for the left sinus. In the ascending aorta, the average 
displacement amplitude is 0.54 mm in the right and noncoronary portion and 
0.34 mm in the left portion of the aortic wall. 

The global distribution of the maximum principal stress was similar in the 
symmetric and asymmetric models, with the maximum stress being concen¬ 
trated at the sinotubular junction, in correspondence with the sinuses; see 
fig. 15. Differences in the maximum principal stress at the same locations for 
the right, left, and noncoronary sinus, are more complex. Whereas the maxi¬ 
mum principal stress in the aorta is essentially independent of the sector of the 
geometry considered, there are large differences in the stress computed in the 
left, right, and noncoronary sinus, especially in the cusp of the sinuses and in 
the sinotubular junction. Peak maximum principal stress for the sinotubular 
junction is 180 kPa, a value comparable to the 220 kPa value obtained by 
Conti et al. [15]. Peak principal stresses in the sinuses, evaluated at the cusp, 
resulted smaller in magnitude than the values reported by Grande-Allen et 
al. [33] , although the relative differences between the stresses in the sinuses is 
similar to those observed in that work. 
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Fig. 13: Results of a dynamic fluid-structure interaction analysis of the asym¬ 
metric aortic root over the final three simulated cardiac cycles, (a) Prescribed 
left-ventricular driving pressure and computed aortic loading pressure, and (b) 
computed flow rate, in which the dashed vertical lines highlight the moment 
that the valve leaflets open and the solid vertical lines highlight the moment 
that the valve leaflets close. Large oscillations indicate the reverberations of 
the aortic valve leaflets upon closure. 
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Fig. 14: Valve leaflet configurations during (a) valve opening and (b) shown 
at equally spaced times during the fourth simulated cardiac cycle. Panels (c) 
and (d) show the corresponding axial velocity profiles. Panel (e) shows the 
computed flow rates, with the dashed lines indicating the times shown during 
valve opening and with the solid lines indicating the times shown during valve 
closure. 


4 Discussion 

In this work, previous IB models of aortic valve dynamics [35,39] are substan¬ 
tially extended by incorporating a deformable model of the aortic sinuses and 
ascending aorta. This extended model is thereby able to more closely replicate 
the complex in vivo dynamics of the real aortic root. As before, the dynamic 
analyses of the aortic valve included four complete cardiac cycles, including a 
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Fig. 15: Displacement contours along the axial section of the asymmetric aortic 
root during valve opening and closure, shown at the same time instants as in 
fig. 14, (a) during valve opening and (a) during valve closure. 


single initialization cycle, which are shown to be sufficient to yield essentially 
periodic results. As in previous IB models of aortic valve dynamics [35,39], 
key hemodynamic parameters, including aortic pressure, maximum flow rate, 
cardiac output, and stroke volume, are in reasonable agreement with physio¬ 
logical ranges [23,77] and are also in good agreement with the patient-specific 
data [56] used to fit the Windkessel model [70] employed in this study. Further, 
grid convergence results demonstrate for the first time that the IB method is 
able to yield essentially grid-resolved bulk flow parameters and vessel defor¬ 
mations at practical grid spacings. 

Many strategies have been developed to determine the zero-pressure ge¬ 
ometry of blood vessels [1,6,21,31,66,74]. Indeed, medical images of arteries 
always provide the vessel geometry in a loaded configuration. It is known, 
however, that using the imaged geometry as the reference geometry of the 
vessel underestimates the amount of strain present in the blood vessel [21]. 
Among the approaches developed to estimate the zero-pressure geometry of 
blood vessels, we chose an iterative backward displacement approach, as de¬ 
scribed by Bols et al. [6] and by Sellier [66]. Iterative approaches based on 
geometric considerations are straightforward to implement and can be applied 
to a variety of problems and solvers because they require minimal implementa¬ 
tion effort. In this work, the backward displacement analysis was implemented 
in ABAQUS and converged to a zero-pressure geometry that, when loaded, 
yielded a deformed configuration that is in good agreement with the initially 







34 


Vittoria Flamini et al. 






(a) (b) 




(c) 


(d) 


Fig. 16: Displacements of the asymmetric aortic root model at four locations: 
(a) the straight portion of the aorta; (b) just above one of the commissures; 
(c) the sinotubular junction; and (d) the middle of the sinus. 


prescribed loaded geometry; see figs. 3-5. However, the deflated geometry is 
quite different from the configuration assumed by a deflated aortic root har¬ 
vested from a body; see fig. 3. This is a drawback of these methods, as shown 
also in a recent study by Wittek et al. [78], in which the unloaded configura¬ 
tion of the aorta is different from an unloaded aorta harvested from a body. 
This mismatch could be the result of many factors, such as the choice of the 
boundary conditions or the need for adding the residual stress to balance the 
level of stress within the aortic wall [28] . We expect that the addition of resid¬ 
ual stress estimates, planned in future work, will reduce the mismatch between 
the deflated geometry computed and the deflated aortic root observed in vivo. 

Interestingly, the addition of a compliant aortic wall did not substantially 
alter the global hemodynamic parameters previously determined by a non- 
compliant IB aortic root model [35]. In general, we observed a slight increase 
in the maximum flow rate as compared to the noncompliant model, whereas 
changes in other global flow parameters were relatively minor. This finding 
might suggest that our vessel wall model is overly stiff despite our use of a 
constitutive model fit to biaxial tensile test data from human aortic tissues 
samples. In contrast, we observed that using a compliant aorta affected more 
clearly the dynamics of valve leaflet closure. Figs. 7(b) and 8(b) show the flow 
inversion and the distension in the sinus, which collects fluid during systolic 

































Immersed boundary-finite element model of FSI in the aortic root 


35 




Max Principal Stress (kPa) 



Fig. 17: Maximum principal stress contours along the axial section of the 
asymmetric aortic root during valve opening and closure, shown at the same 
time instants as in fig. 14, (a) during valve opening and (b) during valve 
closure. 




ejection, and the elastic recoil of the sinus, which pushes down on the leaflet 
during closure. The addition of this mechanism made the model substantially 
more robust. Specifically, we found that the rigid aortic root model required 
the valve leaflet properties to be very carefully tuned in order for the model 
to remain competent throughout the cardiac cycle. This is in clear contrast to 
the present compliant aortic root model, which yields a competent valve for a 
relatively broad range of leaflet parameters. 

Similar model predictions were obtained for both Cartesian grid spacings 
used in this study. We observed hemodynamic values that were within 5% of 
the clinical values reported by Murgo et al. [56] except for the peak flow rate, 
which was within 5% of the values reported by Murgo et al. on the coarser 
Cartesian grid and within 15% on the finer grid. Cardiac output and stroke 
volume were both in very good agreement with the clinical data, within 4% on 
the coarser Cartesian grid, and within 1% on the finer grid. In addition, the 
pulsatile radius, a quantity used to measure aortic stiffness, was also shown to 
be in relatively good agreement with literature values in both cases [42]. 

To make the description of the aortic root dynamics more realistic, we 
modified the aortic geometry and the leaflets configuration to replicate the 
physiological asymmetry of the sinuses and the leaflets. Previous work by 
Grande-Allen et al. [33] and by Conti et al. [15] showed that the use of a realis¬ 
tic geometry in FE and computational fluid dynamics simulations of the aortic 
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Fig. 18: Maximum principal stress of the asymmetric aorta over the cycles that 
reached periodic steady state at four locations: (a) the straight portion of the 
aorta; (b) just above the commissures; (c) the sinotubular junction; and (d) 
the middle of sinus. 


root yields more accurate results for the stress distribution along the sinuses. 
We based our asymmetric aortic root geometries on the inter-commissures dis¬ 
tances reported by Berdajs [5]. Overall, the bulk flow parameters were closer 
to the original values reported by Murgo [56] than the corresponding values 
obtained with a symmetric model for the same level of grid refinement. These 
results suggest that the use of a more realistic geometry could yield more 
accurate estimates of key hemodynamic quantities. Structural stresses gener¬ 
ated by the immersed boundary model were also in good agreement with the 
stresses reported by Conti et al. [15]. Moreover, the level of stress and the 
displacement of the right and in the noncoronary sectors of the aortic root, 
but not in the left sector of the aortic root, were similar in the symmetric and 
the asymmetric case; see figs. 10 and 16. 

Leaflet dynamics for the symmetric aortic valve are shown in figs. 7 and 
8, and leaflet dynamics for the asymmetric aortic valve are shown instead in 
fig. 14. It is possible to see, however, that in the asymmetric model, the left 
coronary leaflet starts closing before the other two leaflets, thus bending the 
flow profile in the aortic tract; see fig. 14(d). Because the aortic root model is 
asymmetric, there is no reason to expect that the dynamics of the three leaflets 
will be identical. We plan to study this asymmetry in the leaflet kinematics 
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more carefully in the context of a more geometrically realistic model of the 
aortic root that accounts for the curvature of the real aortic root and ascending 
aorta. 

4.1 Limitations and Future Work 

Work is underway to address some of the limitations of these models. As de¬ 
scribed previously, our present models do not account for the curved geometry 
of the aortic root, or for the residual stresses present in real arterial vessels. 
We specifically aim to incorporate a description of residual stresses and realis¬ 
tic medical imaging-derived anatomic geometries into our IB models. Further, 
the present models use an isotropic model of the aortic sinuses and ascending 
aorta. Although the experimental data that were used to fit this model show an 
essentially isotropic material response, it is known that real aortic tissues have 
a well-defined collagen fiber structure with an anisotropic material response. 
We also aim to incorporate a more realistic hyperelastic model of the aortic 
wall that accounts for families of collagen fibers [30] as well as experimentally 
constrained models of the elasticity of the aortic valve leaflets [51]. Further, 
it is now well known that the valve leaflets are in fact layered structures [64] . 
These details are not accounted for in the present model, and we anticipate 
that including them in future work will require an important extension of the 
methods presented in this study. Specifically, because the valve leaflets are 
thin elastic structures, we expect that using realistic hyperelastic continuum 
models to describe the mechanics of the valve leaflets will require adopting a 
nonlinear shell formulation. To our knowledge, continuum shell models have 
not yet been widely used within the framework of the IB method, and we an¬ 
ticipate additional computational research will be required to implement these 
types of formulations into our IB simulation framework. 

Another limitation of this work is the long run times required by these 
simulations. Because of the small time step sizes required by the present algo¬ 
rithm, each cardiac cycle required 50,000 or more time steps, requiring on the 
order of one day per cycle on 32-64 cores of a modern high-performance com¬ 
puting (HPC) cluster for the coarse simulations and on the order of one week 
per cycle for the fine-grid simulations. We are currently working to develop 
effective multigrid solvers for stable implicit time stepping schemes for the 
IB method [41]. These methods promise to allow for substantially larger time 
step sizes as well as shorter run times at current grid resolutions. They could 
also enable higher resolution simulations that promise to resolve the fine-scale 
details of the dynamics of the aortic root. 


5 Conclusions 

This work has presented new fluid-structure interaction models of the aortic 
root that extend earlier IB models of the aortic valve by incorporating a finite- 
strain model of the aortic root that uses a constitutive model fit to tensile test 
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data obtained from human aortic root tissue specimens. The models capture 
both the complex fluid dynamics of the flow and the finite deformations of the 
vessel wall, and they yield results that are in good agreement with the clinical 
data that were used to fit the circulation model used in this study, as well 
as to physiological data in the research literature. In addition, a grid conver¬ 
gence study demonstrates that nearly grid converged results are obtained at 
practical spatial resolutions, although yet higher spatial resolution is needed 
to obtain fully resolved simulation results. Finally, by comparing the dynamics 
of symmetric and asymmetric aortic root models, we find that the asymmet¬ 
ric model appears to yield a more accurate description of both the bulk flow 
properties and also the aortic wall mechanics. 
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